Pharmacokinetic Simulation and Area under the Curve Estimation of Drugs Subject to Enterohepatic Circulation

Enterohepatic circulation (EHC) is a complex process where drugs undergo secretion and reabsorption from the intestinal lumen multiple times, resulting in pharmacokinetic profiles with multiple peaks. The impact of EHC on area under the curve (AUC) has been a topic of extensive debate, questioning the suitability of conventional AUC estimation methods. Moreover, a universal model for accurately estimating AUC in EHC scenarios is lacking. To address this gap, we conducted a simulation study evaluating five empirical models under various sampling strategies to assess their performance in AUC estimation. Our results identify the most suitable model for EHC scenarios and underscore the critical role of meal-based sampling strategies in accurate AUC estimation. Additionally, we demonstrate that while the trapezoidal method performs comparably to other models with a large number of samples, alternative models are essential when sample numbers are limited. These findings not only illuminate how EHC influences AUC but also pave the way for the application of empirical models in real-world drug studies.


Introduction
Enterohepatic circulation (EHC) represents a dynamic kinetic process integral to the pharmacokinetics of certain drugs, delineating a crucial pathway within the body's metabolic landscape [1].This intricate cycle commences with drug absorption from the intestine, where a portion of the dissolved drug traverses through the portal system to reach the hepatic lobules.Here, within the intricate architecture of the liver, the drug navigates through branches of the hepatic artery and portal vein, eventually circulating to the central vein.Along this journey, a fraction of the drug diffuses through the sinusoidal walls, finding its way into hepatocytes [2,3].Once within these hepatic cells, the drug may undergo active transportation, facilitated by an array of pumps, such as BCRP, MATE1, MDR3, or MRP2, towards the bile canaliculi, ultimately contributing to bile formation [2,3].Subsequently, bile, laden with the drug, travels through the bile ducts to the gallbladder, where it awaits the signal for release prompted by food intake and the hormone cholecystokinin [4,5].The gallbladder then empties its contents, a process spanning approximately 60 min, through the sphincter of Oddi into the duodenum [6][7][8].Within the duodenum, the drug dissolved in bile may undergo reabsorption, initiating a new cycle of circulation.Figure 1a elucidates these intricate processes alongside parallel kinetic phenomena.
In studying EHC, researchers analyze time versus plasma drug concentration curves to uncover the dynamic interplay between absorption, distribution, metabolism, and excretion, revealing peaks and troughs corresponding to distinct phases of the enterohepatic cycle.Despite offering valuable insights, such analysis is challenged by variability in physiology, drug properties, and environmental factors, along with practical considerations like sampling frequency and timing.Understanding EHC is crucial in drug studies, as it profoundly influences pharmacokinetic parameters essential for determining pharmaceutical efficacy and safety.Consequently, accurately characterizing these plasma profile fluctuations poses a significant challenge in EHC studies, necessitating the development of precise models.The yellow circles represent the drug transfer parameters between the compartments, the majority being first-order constants.Green ellipses represent drug transfer parameters that are defined by an equation.
In studying EHC, researchers analyze time versus plasma drug concentration curves to uncover the dynamic interplay between absorption, distribution, metabolism, and excretion, revealing peaks and troughs corresponding to distinct phases of the enterohepatic cycle.Despite offering valuable insights, such analysis is challenged by variability in physiology, drug properties, and environmental factors, along with practical considerations like sampling frequency and timing.Understanding EHC is crucial in drug studies, as it profoundly influences pharmacokinetic parameters essential for determining pharmaceutical efficacy and safety.Consequently, accurately characterizing these plasma profile fluctuations poses a significant challenge in EHC studies, necessitating the development of precise models.
Physiologically based pharmacokinetic models represent the pinnacle of complexity within EHC modeling, treating organs and tissues as distinct entities with specific properties for both organism and drug [30][31][32][33].However, despite their increasing prevalence, their translation into clinical therapeutics remains limited [34,35].
Alternatively, empirical mathematical functions, devoid of physiological or compartmental underpinnings, have been employed in EHC modeling [36].The pioneering model in EHC employed two gamma functions to fit the plasma mycophenolic acid profile, while another model utilized the sum of two Gaussian functions to depict the plasma profile of simvastatin.Demonstrated feasibility extends these models beyond two peaks by incorporating additional terms [37][38][39].While these models may not elucidate EHC's underlying processes, they align with experimental data and capture the secondary peak characteristics of EHC.Consequently, they have found widespread implementation in clinical settings as valuable tools in therapeutic drug monitoring (TDM).
Understanding the magnitude of EHC for drugs assumes paramount importance, as it significantly influences key pharmacokinetic parameters such as clearance (Cl), volume of distribution (Vd), half-life plasma (t 1/2 ), and AUC [40].Precisely, the AUC is frequently used in TDM for dose adjustment in drugs with a narrow therapeutic range [41].Sampling times and the number of samples are a limiting factor in TDM that requires restricting the sampling strategy.On the other hand, the AUC is an essential parameter in bioequivalence studies, and it has been argued that the conventional trapezoidal method, used to calculate the AUC, probably overestimates its value when EHC exists.
The debate over whether empirical models offer advantages over traditional methods, like the trapezoidal rule, in estimating AUC remains unresolved.Thus, this study aims to evaluate the efficacy of several empirical models in characterizing time-concentration profiles of hypothetical and real drugs sharing the trait of undergoing EHC.This evaluation, grounded in simulations of various pharmacokinetic scenarios, aims to shed light on the potential of empirical models in understanding and quantifying EHC's impact from both conventional and limited sampling perspectives.

Materials and Methods
A compartmental model and a set of empirical models were developed.The compartmental model was used to perform simulations of different pharmacokinetic scenarios; the empirical models were used to estimate the AUC of these simulations.

Compartmental Pharmacokinetic Model
The compartmental model was built considering the kinetic processes involved in EHC, in addition to those that can occur in parallel and that contribute to global pharmacokinetics; see Figure 1a.
The behavior of the EHC phenomenon is schematized in Figure 1b and is given the following system of differential equations: where D x is the concentration of drug and whose subscript indicates the site in which it is located (D c central compartment; D P peripheral; D G intestinal; D B biliary; D S in the solid pharmaceutical form).Except for the parameters f dis , k gem , and v bio , which are set as parameters as a function of time or drug concentration, most of the parameters in this system of equations are constant during the simulation.
Fraction of dissolved dose of solid dosage form ( f dis ) is one of the time-dependent parameters and was simulated with the Weibull function in Equation ( 2).This is a function without a physicochemical basis widely used in modeling drug dissolution [42], and it is expressed as follows: where t lag corresponds to the time it takes to start the dissolution process, k d is the scale parameter that determines the fraction of drug released from the formulation, and a is the shape parameter that characterizes the curve as exponential (a = 1), sigmoid (a > 1), or parabolic (a < 1).
Emptying of the gallbladder (k gem ) is a time-dependent parameter; continuous function, given by Equation (3), was used to simulate emptying of the gallbladder and allows the generation of several emptying cycles distributed over a period of 24 h; it is expressed in the following way: where the parameter k max is the maximum emptying speed of the gallbladder, the parameter b is related to the duration of each of the emptying cycles, and the parameter c i allows the establishment of the moments in which the emptying of gallbladder occurs for each of the EHC cycles according to meal times (established as hours elapsed after drug intake).
The last non-constant parameter v bio is a function of the drug concentration in the central compartment and represents the biotransformation rate of the drug.For its calculation, the Michaelis-Menten equation, Equation (4), is used.This equation allows us to consider the effects of the competition that may exist between biotransformation and other processes that the drug may have (bile transfer, transfer to the peripheral compartment, or renal elimination), and it is expressed as follows: where k m is the Michaelis-Menten constant, and v max is the maximum biotransformation rate.
To complete the analysis, the percentage of drugs within the central compartment undergoing enterohepatic circulation was calculated.This percentage is not necessarily constant over time; Equation (5) allows the calculation of its temporal fluctuation, while Equation (6) calculates its theoretical maximum value. (5)

Compartmental Model Assumptions
The compartmental model was developed with the objective of representing the EHC and some parallel pharmacokinetic processes in the most realistic way possible; however, to simplify the analysis, some reasonable assumptions are made: i. Compartments are homogeneous, so drug concentration within each of the compartments reaches an instantaneous kinetic equilibrium; ii.Liver is located within the central compartment, so there is an instantaneous kinetic equilibrium between plasma and intrahepatic drug concentrations; iii.Elimination by routes other than renal, fecal, or biotransformation is insignificant; iv.Fraction of drug that is biotransformed does not need to be reconverted to the original drug to undergo EHC, as this usually happens under the presence of microbiota glucuronidases; v. Elimination by biotransformation occurs through a single pathway involving a single enzyme; vi.Drug dissolved in the bile is excreted to the intestinal compartment only from the bile compartment, so direct excretion from the central compartment is non-existent; vii.Fraction of drug dissolved, transfer of drug from the gallbladder to the intestine, and rate of biotransformation can be realistically represented by equations, while the rest of the parameters follow first-order kinetics that do not consider changes in transport rates between compartments due to physiological phenomena.
Finally, during simulations, it was assumed that the gallbladder is only activated by meals three times during a 24-hour period; at other times, its emptying is null.

Pharmacokinetic Simulations
The compartmental model was implemented through the SimBiology tool of MAT-LAB version R2021a.Drug concentrations versus time were simulated with the sundials solvertype.The description, units, and specific values of model parameters are shown in Table 1.These values were proposed to have a starting model for subsequent scenarios.
Data shown in Table 1 correspond to those of the standard model from which the values of one parameter were varied at a time to generate 18 different scenarios.These scenarios are grouped into pairs that differ in whether the value of the modified parameter is greater or less than its value in the standard model.Table 2 shows the description, the modified parameter, and its value in each scenario.To increase the diversity of the data generated in each scenario, the values of the parameters k m and k ehc , which are related to the magnitude of the EHC, were varied.Table 3 shows the values of k m and k ehc for each of the cases and their effect on the EHC.Different combinations of these two parameters allowed the generation of up to nine different simulations per scenario, giving a total of one hundred and sixty-two simulations.At the end of each simulation, post-simulation calculations were performed to obtain the AUC 0−∞ .From each simulation, time versus plasma concentration data were extracted at specific times to represent the sampling process.The sampling schemes were divided into three categories, which are shown in Table 4.

Sampling Type Sampling Times (h) Number of Samples Description
Meal-based sampling 0, 0.
Focuses on sampling before and after food intake Conventional sampling 0, 0.25, 0.

Empirical Pharmacokinetic Models
The empirical models were fitted to the data extracted from each simulation.Five models (M1 to M5) represented by functions were built to characterize the multiple peaks present in the EHC.They were compared with each other and with respect to the reference model of extravascular administration (M0), Equation (7).Equations ( 8) and ( 9) are modifications of previously reported models [37,38], while the rest, Equations ( 10)- (12), are new models built intuitively.They are described as: where the letters a, b, c, d, etc., are the parameters of the models; t is the time and C M (t) the plasma concentration estimated by each model.Except for the model given by Equation ( 1), the models are built with summations that describe the first absorption peak (when i = 1) and secondary resorption peaks caused by EHC (when i > 1); therefore, in these models, the number of peaks is equal to n.Another characteristic that all models have in common is that they incorporate a term or factor to the left of the sum and that characterizes the elimination process.
The model described by Equation ( 12) is an extension of the model in Equation ( 11) and is intended to predict the presence of additional peaks that may not be detected after the first 24 h.This prediction is made under the assumption that the gallbladder will empty with an intensity and at a time like the first day.In this model, the intensity and moment of gallbladder emptying are given by the parameters b i and d i , respectively.The prediction after 24 h is based on the parameters that describe the largest peak (usually the first absorption peak, b i and d i ) and the last observed peak (b n and d n ).
The AUC 0−∞ of all models converge to a finite value when the model is fitted to time versus concentration data; its value can be calculated exactly using explicit formulas.Equations (S1)-(S6) shown in the Supplementary Materials (Supplementary Data S1) allow us to perform this calculation exactly, and they were obtained with the help of Maple software version 2018.

Compilation of Pharmacokinetic Profiles from the Literature
A bibliographic search and compilation of pharmacokinetic profiles of drugs that undergo EHC reported in different works was carried out.The search criteria were drugs administered orally, in a single dose, and with a profile with at least two clearly visible peaks.WebPlotDigitizer version 4.5 software was used to recover time and plasma concentration data from the graph reported in the study.The software allows you to select the points in the pharmacokinetic profile to estimate the values of the time-concentration data pairs.This software has been used successfully for data collection from graphs [43].

Data Fitting, Statistical Analysis, and Evaluation of Empirical Models
Non-linear regression (NLR) of the simulated data retrieved from the literature was performed using the Curve Fitting tool of MATLAB version R2021a.When necessary, restrictions were applied to the coefficients to obtain a good parametric fit.The adjusted R 2 , SSE, and RMSE error metrics were obtained directly from said tool.The script used to execute the NLR is provided in the Supplementary Materials (Supplementary Data S2).
Evaluation of models was carried out with the corrected Akaike Criterion (AICc) for small sample sizes using Equation ( 13): Pharmaceutics 2024, 16, 1044 where k is the number of model parameters, and n is the number of observations.As part of evaluation of the models, comparisons were made between the reference AUC 0−∞ of the simulation, with the AUC 0−∞ calculated by the trapezoidal method and estimated by empirical models, Equations (S1)-(S6).

Effect of EHC Degree on AUC
Initially, we aimed to verify whether the compartmental model could accurately simulate the EHC process, allowing us to investigate its impact on AUC.Simulations following single-dose extravascular administration delineated the relationship between gallbladder emptying, gallbladder drug concentrations, and plasma profile (Figure 2).Notably, spikes in plasma profile coincided with gallbladder contractions, underscoring the dynamic nature of EHC.Additionally, drugs subject to extensive EHC (depicted by the black line) exhibited prolonged residence in the organism compared to those experiencing minimal EHC (illustrated by the yellow line).Expanding our model to encompass biotransformation and gallbladder transport enabled us to assess the influence of their interplay.Simulations with a  = 50 revealed a consistent EHC% value (bold lines, left panel of Figure 3), suggesting negligible competition.Conversely, simulations with  = 0.005 yielded fluctuating EHC% values, ranging between a minimum and a theoretical maximum determined by Equation (6) (pale lines, right panel of Figure 3), indicative of significant competition.In this scenario, maximum EHC% values ranged from 88% to 33%.This competitive dynamic directly impacted  , where cases of negligible competition yielded higher maximum values (ranging from 20.2 to 30.9 mg × h/L) compared to scenarios with relevant competition, where  decreased with decreasing  values.These findings affirm the completeness and realistic simulation capabilities of our model to simulate the EHC process.Expanding our model to encompass biotransformation and gallbladder transport enabled us to assess the influence of their interplay.Simulations with a k m = 50 revealed a consistent EHC% value (bold lines, left panel of Figure 3), suggesting negligible competition.Conversely, simulations with k m = 0.005 yielded fluctuating EHC% values, ranging between a minimum and a theoretical maximum determined by Equation (6) (pale lines, right panel of Figure 3), indicative of significant competition.In this scenario, maximum EHC% values ranged from 88% to 33%.This competitive dynamic directly impacted AUC, where cases of negligible competition yielded higher maximum values (ranging from 20.2 to 30.9 mg × h/L) compared to scenarios with relevant competition, where AUC decreased with decreasing k m values.These findings affirm the completeness and realistic simulation capabilities of our model to simulate the EHC process.

Simulation of Pharmacokinetic Scenarios and Fitting to Empirical Models
We conducted 162 simulations across eighteen scenarios to explore their impact on  , as summarized in Table 5. Comparative analysis revealed that scenarios S1 through  and  exhibited significantly different  compared to standard model  .Additional details and pharmacokinetic simulation plots for each scenario are available in the Supplementary Materials (Supplementary Data S3, Figures S1-S18).Next, we fitted six empirical models to the obtained plasma profile data.Figure 4 illustrates the fit of these models to the data from a representative simulation.Apart from model 0, all models successfully captured the multiple peaks of the plasma profile, exhibiting subtle visual differences.

Simulation of Pharmacokinetic Scenarios and Fitting to Empirical Models
We conducted 162 simulations across eighteen scenarios to explore their impact on AUC 0−∞ , as summarized in Table 5. Comparative analysis revealed that scenarios S1 through S 6 and S 15 exhibited significantly different AUC 0−∞ compared to standard model S 7 .Additional details and pharmacokinetic simulation plots for each scenario are available in the Supplementary Materials (Supplementary Data S3, Figures S1-S18).Next, we fitted six empirical models to the obtained plasma profile data.Figure 4 illustrates the fit of these models to the data from a representative simulation.Apart from model M0, all models successfully captured the multiple peaks of the plasma profile, exhibiting subtle visual differences.
Criteria for selecting the best model included adjusted R 2 , AICc, and accuracy in predicting AUC. Figure 5 displays scatter plots of these criteria, while Table 6 presents their means for each model.Figure 5 allows the visual identification of the result of each regression of the empirical models in a scatter plot.Models M4 (cyan balls) and M5 (red balls) stand out for showing a higher degree of clustering than other models, and both models have the lowest AICc values and the highest adjusted R 2 values.The gamma distribution model (pink balls) presented poorer clustering due to the greater difficulty of performing reproducible regressions since, in this model, applying restrictions for a successful NLR was less intuitive for the authors.Criteria for selecting the best model included adjusted  , AICc, and accuracy in predicting .Figure 5 displays scatter plots of these criteria, while Table 6 presents their means for each model.Figure 5 allows the visual identification of the result of each regression of the empirical models in a scatter plot.Models 4 (cyan balls) and 5 (red balls) stand out for showing a higher degree of clustering than other models, and both models have the lowest AICc values and the highest adjusted  values.The gamma distribution model (pink balls) presented poorer clustering due to the greater difficulty of performing reproducible regressions since, in this model, applying restrictions for a successful NLR was less intuitive for the authors.
Model 5 achieved the highest values for adjusted  (0.979),  (−58.1), and accuracy in predicting  (105.9).Although 5 exhibited the smallest spread in adjusted  (±0.02), it also had the highest spread in  prediction accuracy (±21.4).These results indicate that 5 demonstrated the best performance characteristics and was thus chosen for subsequent analysis.predicting .Figure 5 displays scatter plots of these criteria, while Table 6 presents their means for each model.Figure 5 allows the visual identification of the result of each regression of the empirical models in a scatter plot.Models 4 (cyan balls) and 5 (red balls) stand out for showing a higher degree of clustering than other models, and both models have the lowest AICc values and the highest adjusted  values.The gamma distribution model (pink balls) presented poorer clustering due to the greater difficulty of performing reproducible regressions since, in this model, applying restrictions for a successful NLR was less intuitive for the authors.Model 5 achieved the highest values for adjusted  (0.979),  (−58.1), and accuracy in predicting  (105.9).Although 5 exhibited the smallest spread in adjusted  (±0.02), it also had the highest spread in  prediction accuracy (±21.4).These results indicate that 5 demonstrated the best performance characteristics and was thus chosen for subsequent analysis.Model M5 achieved the highest values for adjusted R 2 (0.979), AICc (−58.1), and accuracy in predicting AUC (105.9).Although M5 exhibited the smallest spread in adjusted R 2 (±0.02), it also had the highest spread in AUC prediction accuracy (±21.4).These results indicate that M5 demonstrated the best performance characteristics and was thus chosen for subsequent analysis.

Estimation of AUC in Different Pharmacokinetic Scenarios
Applicability of model M5 to various pharmacokinetic scenarios was examined.Results indicated that model M5 accurately predicted AUC for most scenarios, with close to 100% accuracy (Figure 6).However, scenarios S 3 (rapid renal elimination from central compartment) and S 15 (fast peripheral distribution with slow central return) notably overestimated the actual AUC (exceeding 110% and 150%, respectively).This overestimation persisted regardless of whether the trapezoidal method or model M5 was employed.

Estimation of AUC in Different Pharmacokinetic Scenarios
Applicability of model 5 to various pharmacokinetic scenarios was examined.Results indicated that model 5 accurately predicted  for most scenarios, with close to 100% accuracy (Figure 6).However, scenarios  (rapid renal elimination from central compartment) and  (fast peripheral distribution with slow central return) notably overestimated the actual  (exceeding 110% and 150%, respectively).This overestimation persisted regardless of whether the trapezoidal method or model 5 was employed.

Estimation of AUC under Different Sampling Schemes
Next, we investigated the impact of sampling schemes on predictions of model 5.Firstly, we contrasted meal-based and conventional sampling methods (Figure 7).Notably, when the percentage of drugs undergoing EHC was 40%, model-predicted  accuracy was nearly 100% for both sampling schemes.However, at 80% EHC, accuracy decreased to approximately 50% with conventional sampling, while remaining close to 100% with meal-based sampling.Regardless of the sampling scheme or the percentage of drugs undergoing EHC,  estimates by the trapezoidal method closely mirrored those by the model.
Secondly, different sampling schemes for TDM were assessed (Figure 8).It was observed that with TDM1 and TDM2 schemes, model-predicted  accuracy ranged between 50-75%, depending on the percentage of drugs undergoing EHC.In contrast, TDM3 and TDM4 schemes notably improved accuracy to over 75%, while TDM5 showed a modest increase.Notably,  predictions by the model were nearly identical with TDM1/TDM2 and TDM3/TDM4 schemes, respectively.In all cases,  estimates by the trapezoidal method improved gradually with increased sampling times.

Estimation of AUC under Different Sampling Schemes
Next, we investigated the impact of sampling schemes on predictions of model M5.Firstly, we contrasted meal-based and conventional sampling methods (Figure 7).Notably, when the percentage of drugs undergoing EHC was 40%, model-predicted AUC accuracy was nearly 100% for both sampling schemes.However, at 80% EHC, accuracy decreased to approximately 50% with conventional sampling, while remaining close to 100% with mealbased sampling.Regardless of the sampling scheme or the percentage of drugs undergoing EHC, AUC estimates by the trapezoidal method closely mirrored those by the model.
Secondly, different sampling schemes for TDM were assessed (Figure 8).It was observed that with TDM1 and TDM2 schemes, model-predicted AUC accuracy ranged between 50-75%, depending on the percentage of drugs undergoing EHC.In contrast, TDM3 and TDM4 schemes notably improved accuracy to over 75%, while TDM5 showed a modest increase.Notably, AUC predictions by the model were nearly identical with TDM1/TDM2 and TDM3/TDM4 schemes, respectively.In all cases, AUC estimates by the trapezoidal method improved gradually with increased sampling times.4.  4. Figure 8. Fit of model M5 to drug profiles submitted to 40% (top), 60% (middle), and 80% (bottom) of EHC%.Adjustment to the data was performed by sampling for the TDM, as shown in Table 4.
The AUC was estimated using model M5 and compared with reported values, as shown in Figure 9, obtaining a correspondence between both values close to 100%.The largest discrepancies were observed for Ezetimibe and Meloxicam, with differences of 78.4% and 124.1%, respectively, while other estimates fell within the range of 80-120%, as depicted in Table 7.
The  was estimated using model 5 and compared with reported values, as shown in Figure 9, obtaining a correspondence between both values close to 100%.The largest discrepancies were observed for Ezetimibe and Meloxicam, with differences of 78.4% and 124.1%, respectively, while other estimates fell within the range of 80-120%, as depicted in Table 7.  a The elimination half-life of amiodarone is biphasic; the half-life of the first phase of disposition is presented here since the actual half-life is between 25 and 110 days.a The elimination half-life of amiodarone is biphasic; the half-life of the first phase of disposition is presented here since the actual half-life is between 25 and 110 days.

Effect of EHC Degree on AUC
Simulations of our compartmental model are in agreement with other studies in which it is stated that the EHC increases the permanence of the drug in the body by increasing t 1/2 apparent elimination [40,57].We want to highlight the linear nature of our compartmental model in which an increase in dose proportionally increases AUC regardless of the degree of EHC% (Supplementary Data S5, Figure S29).This is not surprising considering that drugs such as atorvastatin and mycophenolic acid that experience EHC maintain linear pharmacokinetics with respect to AUC [58,59].However, a less understood effect of the EHC is its impact on AUC.We take advantage of the information generated from the simulations to address a question that has been discussed in other works.
It has been discussed on several occasions how AUC is altered when the drug undergoes EHC.The prevailing position is that AUC increases when EHC exists; several papers have provided evidence to support this position [60][61][62][63].However, in other studies, it has been proposed that the degree of EHC does not have an impact on AUC [57,64].Recently, Ibarra et al. reported that the inclusion of a hepatic compartment, where the competition between hepatobiliary secretion and biotransformation is explicit, allows us to explain within their model why the level of AUC increases when the degree of EHC is higher [40].Similar to the model proposed by Ibarra, in our model, this competition is considered with the difference that the elimination constant by biotransformation is not of the first order but follows Michaelis-Menten kinetics.
For drugs with a very large k m , the rate of elimination by biotransformation (v bio ) is negligible, but when the k m rate is very small, v bio increases, and the competition between bile secretion and biotransformation becomes important.When this happens, v max is reached more quickly due to the drug's high affinity for the enzyme.As described above, the competition between the two processes causes a fluctuation of EHC% between a theoretical minimum and a theoretical maximum.The maximum can be estimated with Equation ( 6), while Equation ( 14) allows estimating the minimum possible for EHC%, which is reached if the concentration of drug in the compartment is high enough.
Figure 3 shows that AUC of drug is lower at EHC max % equal to 88% compared to those with EHC max % equal to 33%.Within the proposed model, this is explained by the fact that biliary transport behaves as an additional elimination pathway that results in a decrease in AUC 0−∞ by the amount of the drug extracted from the central compartment.
According to our model, we can differentiate two ways that increase AUC by altering the degree of EHC that a drug undergoes.The first is by decreasing k m of biotransforming enzyme by adding a competitor, and the second is by suppressing drug transport into the gallbladder by decreasing the value of k ehc .We translate these two possibilities into a couple of examples.
In the first case, it could happen that a drug that is subjected to EHC and to biotransformation is administered concomitantly with another drug that is also biotransformed by the same enzyme (e.g., CYP3A4).Competition of both substrates for the same enzyme will increase k m of the drug that EHC is experiencing.As shown in Figure 3, high values of k m increase AUC.In fact, simvastatin, a drug that undergoes EHC, and cyclosporine are biotransformed by the same enzyme CYP3A4; when administered concomitantly, the body's exposure to simvastatin increases enough to give rise to the toxic effects of simvastatin [65].
In the second case, a drug could be presented that undergoes EHC in the normal population, but where a mutation in the gene encoding a protein responsible for drug transport to the biliary canniculus (e.g., MRP2) is able to decrease the activity of the transporter in carriers of the mutation.According to our model, this mutation would imply a decrease in EHC max % and an increase in AUC with respect to the normal population.In fact, people who suffer from the so-called Dubin-Johnson syndrome have higher levels of bilirubin in plasma than a normal person because the MRP2 transporter is defective, preventing bilirubin from experiencing enterohepatic circulation to which it is normally subjected, which causes its accumulation within the body [66].
In summary, suppressing biliary transport or favoring biliary transport over biotransformation could increase AUC of drugs from two different perspectives.In one study, both perspectives were investigated and reached conclusions similar to those of our model [67], but more studies are still needed to investigate how much of the AUC increase observed in these cases can be attributed to the modification of the EHC.

Simulation of Pharmacokinetic Scenarios and Fit to Empirical Models
In addition to the biliary transfer rate (k ehc ), we can point out from Table 5 that the parameters related to absorption (k a ), renal elimination (k rel ), intestinal elimination (k gel ), and transport between compartments (k 12 y k 21 ) were the ones that most influenced the modification of AUC 0−∞ with respect to the standard model.Regardless of the effect that these processes may have, all of them have the characteristic pharmacokinetic profile of "sawtooth", a sufficient requirement for their adjustment to empirical models.
A qualitative analysis of Figure 4 reveals that except for the classical model of extravascular administration (M0), most models can capture the multi-peak phenomenon, although with some differences.For example, the model M1 generated symmetric peaks, while the rest of the models generated asymmetric peaks more similar to those in the simulation.The models M1 and M2 presented a smooth transition from one peak to the next, while in the rest of these models, the transition was more abrupt.The models M1 and M2 tended to overestimate the magnitude of some of the peaks.The models M4 and M5 generated nearly identical predictions in the first 24 h, but M5 stood out for the ability to predict subsequent peaks even if they were not identified by the sampling method.
On the other hand, a quantitative analysis of the results in Table 6 shows that increasing the number of model parameters does not necessarily imply an adjusted R 2 value closer to 1. Models that required more parameters (M1, M2, and M3) had lower R 2 values than models that required fewer parameters (M4 and M5).However, the number of parameters increases if you increase the value of AICc.In fact, model M1, which requires up to 14 parameters to fit a profile with four peaks, had the highest AICc value of all the models evaluated.As complex as it may seem, model M5, described by Equation ( 12), requires relatively few parameters to describe various peaks, which is reflected in the smallest AICc of all the models evaluated.
It is important to mention that the greater the number of secondary peaks that you want to model, the greater the number of parameters.Increasing the number of parameters also implies increasing the number of samples; otherwise, the NLR could not be performed.The number of samples must be greater than or equal to the number of parameters to be used to characterize one or more peaks.Table 8 shows the number of parameters (and, therefore, minimum samples required) to characterize the peaks of a plasma profile with EHC.Depending on the number of peaks observed, you might choose to use more or fewer parameters for a given model.Regarding the accuracy of the prediction AUC, it is interesting to note that even the classical model of extravascular administration (M0) that is unable to capture secondary peaks has an accuracy in AUC estimation similar to some of the models that can capture secondary peaks.Previously, we presented an empirical model to capture the secondary peaks of the plasma profile of simvastatin [38], which is equivalent to the model M1 of this work.Unlike our previous work, this time, the biexponential model M0 was better at estimating AUC than the model M1.This suggests that model M1 might fit well with the plasma profiles of simvastatin but not with other drugs with different pharmacokinetic characteristics.
With respect to the model first proposed by Premaud et al. [37] and that in this work, which is equivalent to the model M2, a couple of questions can be said.This model has undoubtedly been the most successful in practice, as it has demonstrated its superiority in the field of TDM over the fixed dose in the monitoring of mycophenonlate mofetil, a drug with a narrow therapeutic margin [68].Adaptations to it have even been proposed for the case of profiles with three peaks, which has allowed its usefulness to be further expanded [39].In this paper, we present the model in its general form to be able to model any number of peaks, Equation ( 9), and we also provide a formula for the exact calculation of AUC 0−∞ from its parameters, Equation (S3).However, in this work, there were models that presented a better performance; this was mainly because other models required a smaller number of parameters to model the same pharmacokinetic profile.In addition, on many occasions, the execution of the NLR was less intuitive for M2 compared to other models such as M3, M4, and M5.

Estimation of AUC in Different Pharmacokinetic Scenarios
Once the superiority of model M5 had been demonstrated, we proceeded to evaluate under which conditions it could be more useful and in which it does not offer significant improvements compared to the trapezoidal method.As mentioned above, the M5 model overestimated AUC in scenarios S 3 and S 15 .This suggests that AUC predictions from model M5 are more likely to be inaccurate for drugs that are eliminated very quickly or those with two-compartment kinetics in which the drug remains for longer in the body.It is important to note that kinetic equilibrium is not reached instantaneously between compartments, and the addition of a peripheral compartment in scenario S 15 makes it possible to take into account the presence of highly lipophilic drugs.Considering this, model M5 might have trouble estimating AUC of highly lipophilic drugs.
An important aspect to highlight is that AUC calculation for each scenario using Equation (S6) was not significantly different from the one calculated by the trapezoidal method.These results suggest that, if possible, the selection of a meal-based sampling scheme is more important than the selection of a complex model that fits the data well but does not offer practical improvements beyond a good fit.
In the past, it has been argued that possibly the use of the trapezoidal method in drugs undergoing EHC could lead to an overestimation of AUC.This assumption is based on the fact that the lack of resolution of the secondary peaks would lead to an overestimation of the actual area of these peaks, similar to approximating the area of a circle inscribed within a square from the area of that square.However, in this paper, we observed that, more often than not, the trapezoidal method estimates the real quite well, at least when there is an exhaustive sampling based on meals; see Figure 6.

Estimation of AUC under Different Sampling Schemes
The results suggest that it is more important to opt for a meal-based sampling method when secondary peaks are very pronounced due to a high degree of EHC than when peaks are very small due to a low degree of EHC.Even so, in practice, it would be worthwhile to opt for meal-based sampling whenever possible, as the contribution of the EHC process to the AUC of some drugs can be highly variable.For example, the variability of EHC in mycophenolic acid ranges from 10 to 60% [69].
It is important to note that in a comprehensive meal-based or conventional sampling, the AUC calculated by the trapezoidal method is very similar to that calculated by model M5; see Figure 7.This result suggests that when sufficiently large numbers of samples are available, as in the case of bioequivalence studies, the differences between each of the methods for calculating AUC are negligible.In these cases, opting for the trapezoidal method would be the most practical.
When sampling was limited, we found that increasing the number of samples does not necessarily lead to an improvement in the estimation of AUC.True improvements in the estimation of AUC were observed between TDM2 to TDM3 and TDM4 to TDM5 schemes.A more in-depth analysis of Figure 8 shows that the TDM1 and TDM2 samples allow the characterization of only two peaks, TDM3 and TDM4 three peaks and TDM5 four peaks of the plasma profile.The clear conclusion is that it is more important to increase the number of properly characterized peaks than to increase the number of samples.
It is important to note that in our study, the sampling times for TDM were not chosen arbitrarily but followed the rationality of meal-based sampling.From Figure 8, it is easy to see that the sampling times coincide almost exactly with the maximums and minimums of the plasma profile.These results suggest that in the TDM of drugs subject to EHC, it is essential to establish a meal schedule for the planning of sampling and that these are taken for a sufficiently long time to characterize several peaks.In this sense, our work supports the suggestions made by Lixuan et al., who studied the effect of mealtimes on AUC of mycophenolic acid [70].

Application of Model M5
Our study has the limitation that it does not have data from individual real-world profiles; we only had access to the average plasma profiles of some drugs published in other studies.Even so, we were able to adjust model M5 to all the profiles collected and obtained, in most cases, the estimates of AUC that were quite close to those reported in their respective articles.In the process, we obtained profiles with dynamics for which we cannot be sure that they approximate the real dynamics since several of these studies did not report or consider the mealtimes for sampling.The adjustment of the proposed model to these profiles is the first step to promote its application in data generated from real-world studies based on rational sampling.
It is clear that one of the deficiencies of model M5 is its inability to estimate the degree of EHC; however, as seen in Table 7, this model estimates the AUC quite well.It is not easy to estimate experimentally the degree of EHC that a drug experiences, and there are few works that report this value.Due to the inability to perform an analysis of the fit of model M5 to the profiles of the drugs studied from these data, we opted to use the elimination half-life t 1/2 to perform this analysis.
We previously pointed out that the M5 model had difficulties in adequately estimating the AUC in those scenarios that considered drugs with very long permanence in the body.In Table 7, we find precisely that model M5 considerably overestimated the AUC of meloxicam and amiodarone, which have a very high t 1/2 (39.5 h and more than 25 days, respectively).At the other extreme, model M5 underestimated the AUC of ezetimibe although it does not have a particularly high t 1/2 (between 13.8 and 15.1 h).
It should be clarified that although the empirical models were applied to data generated from simulations of the EHC process, in theory, they could be applied to any other profile that has multiple peaks caused by mechanisms other than the EHC.This would be possible since, being of an empirical nature, they do not need the determination of kinetic parameters of physiological or physicochemical basis.In this way, they could be applied to plasma profiles that show multiple peaks, such as those caused by pulsatile delivery systems [71,72] (examples of these commercially available pulsatile systems are Pulsincap ® , Ritalin ® , Pulsys ® ) or those caused by physicochemical differences between different regions of the gastrointestinal tract.It has been suggested that differences in the pH of the gastrointestinal tract that ionize mycophenolic acid and not EHC could be responsible for the double-peak phenomenon observed in this drug [73].However, it would be inapplicable to vascular delivery drug profiles that exhibit multiple spikes regardless of the mechanism behind these spikes (e.g., fentanyl, sufentanil, propofol, and phenoperidine [74][75][76][77]).
We highlight the feasibility of applying model M5 as a structural model in population pharmacokinetic models.Many structural models employed in population pharmacokinetic approaches have used physiological or mechanism-based concepts to describe the EHC [78].An example of the success of this type of model contains forty-five structural parameters and simultaneously describes the extent and time course of EHC in rats, dogs, and humans for fimasartan [79].However, empirical models with fewer parameters might present advantages as structural models.The most recent example is the non-linear mixed effect (NLME) approach based on stochastic approximation expectation-maximization for the double gamma distribution model [80].
Our model is unable to predict the effect on the pharmacokinetic profile of variables such as age, gender, body weight, pathological states, etc.In this sense, the inclusion of covariates in our model for the development of further population pharmacokinetic models would allow us to limit its application to specific clinical environments in the same way that has been done for the model proposed by Premaud et al. [37].

Conclusions
In this study, we conducted simulations to assess the impact of EHC on AUC in various scenarios.Our findings indicate that modulating biliary transport can effectively increase AUC.We evaluated five empirical models for characterizing EHC peaks, identifying model M5 as the most efficient due to its simplicity, accurate AUC estimates, and ability to predict peaks beyond sampling times.While model M5 performed well in most scenarios, it struggled with drugs exhibiting rapid renal elimination or bicompartmental kinetics.
Evaluation of different sampling schemes revealed that our proposed model did not consistently outperform the trapezoidal method, particularly with extensive sampling.However, limited sampling for TDM benefited from our model's ability to characterize peaks, yielding superior AUC estimates compared to the trapezoidal method.In conclusion, this study underscored the importance of implementing meal-based sampling schemes and laid the groundwork for applying empirical models in real-world studies of drugs undergoing EHC.

Figure 1 .
Figure 1.(a) Kinetic processes present in the EHC.Additional kinetic processes may occur, such as fecal elimination of a fraction of the drug secreted in the bile or not absorbed, renal elimination, distribution to peripheral tissues, hepatic biotransformation, and release of the pharmaceutical form for drugs for extravascular administration; (b) scheme of the compartmental model implemented for the simulation of the EHC.Rectangles with rounded corners represent the model compartments.The yellow circles represent the drug transfer parameters between the compartments, the majority being first-order constants.Green ellipses represent drug transfer parameters that are defined by an equation.

Figure 1 .
Figure 1.(a) Kinetic processes present in the EHC.Additional kinetic processes may occur, such as fecal elimination of a fraction of the drug secreted in the bile or not absorbed, renal elimination, distribution to peripheral tissues, hepatic biotransformation, and release of the pharmaceutical form for drugs for extravascular administration; (b) scheme of the compartmental model implemented for the simulation of the EHC.Rectangles with rounded corners represent the model compartments.The yellow circles represent the drug transfer parameters between the compartments, the majority being first-order constants.Green ellipses represent drug transfer parameters that are defined by an equation.

Figure 2 .
Figure 2. Pharmacokinetic simulations at different percentages of drugs undergoing enterohepatic circulation.The plasma profile in the middle of the drug concentrations in the gallbladder and below the kinetics of gallbladder emptying are shown above, expressed as variation in the kgem according to Equation (3).

Figure 2 .
Figure 2. Pharmacokinetic simulations at different percentages of drugs undergoing enterohepatic circulation.The plasma profile in the middle of the drug concentrations in the gallbladder and below the kinetics of gallbladder emptying are shown above, expressed as variation in the k gem according to Equation (3).

Figure 3 .
Figure 3.Effect of EHC% on cumulative AUC.On the (left), the temporal dynamics of the percentage of drugs subjected to EHC% is shown, and on the (right), the temporal increase in the cumulative AUC is shown.

Figure 3 .
Figure 3.Effect of EHC% on cumulative AUC.On the (left), the temporal dynamics of the percentage of drugs subjected to EHC% is shown, and on the (right), the temporal increase in the cumulative AUC is shown.

Pharmaceutics 2024 , 23 Figure 4 .
Figure 4. Fitting of empirical models to data from a representative simulated pharmacokinetic profile.In this figure, the adjustment was made to the data obtained through a meal-based sampling, where up to 80% of the drug undergoes EHC.

Figure 5 .
Figure 5. Scatter plots of the accuracy of AUC versus AICc and R 2 of each model (M0-M5).

Figure 4 .
Figure 4. Fitting of empirical models to data from a representative simulated pharmacokinetic profile.In this figure, the adjustment was made to the data obtained through a meal-based sampling, where up to 80% of the drug undergoes EHC.

Figure 5 .
Figure 5. Scatter plots of the accuracy of AUC versus AICc and R 2 of each model (M0-M5).

Figure 5 .
Figure 5. Scatter plots of the accuracy of AUC versus AICc and R 2 of each model (M0-M5).

Figure 6 .
Figure 6.Accuracy of AUC predicted by model M5 and by the trapezoidal method for each pharmacokinetic scenario.Calculation was made from data obtained from a meal-based sampling.

Figure 6 .
Figure 6.Accuracy of AUC predicted by model M5 and by the trapezoidal method for each pharmacokinetic scenario.Calculation was made from data obtained from a meal-based sampling.

Figure 8 .
Figure 8. Fit of model M5 to drug profiles submitted to 40% (top), 60% (middle), and 80% (bottom) of EHC%.Adjustment to the data was performed by sampling for the TDM, as shown in Table4.

Figure 8 .
Figure 8. Fit of model M5 to drug profiles submitted to 40% (top), 60% (middle), and 80% (bottom) of EHC%.Adjustment to the data was performed by sampling for the TDM, as shown in Table4.Figure8.Fit of model M5 to drug profiles submitted to 40% (top), 60% (middle), and 80% (bottom) of EHC%.Adjustment to the data was performed by sampling for the TDM, as shown in Table4.

Figure 9 .
Figure 9.Comparison between reported for drugs that undergo EHC and AUC predicted by model M5.

Figure 9 .
Figure 9.Comparison between reported for drugs that undergo EHC and AUC predicted by model M5.

Table 1 .
Parameters of the standard model.

Table 2 . Modified parameters in pharmacokinetic scenarios. Scenarios Parameter Value k m a k ehc a Description
a Values of k m and k ehc are grouped in different cases represented by Roman numerals and capital letters, respectively; see Table3.

Table 3 .
Values of the constants k m and k ehc and their effect on the percentage of drugs that undergo EHC.
a The effect is observed as a decrease with respect to the maximum value: insignificant decrease: +; moderate decrease: ++; outstanding decrease: +++.
a Due to the greater persistence of drugs in the organism, simulation time was longer than 96 h in these scenarios.b Scenarios with p-values less than 0.05.
a Due to the greater persistence of drugs in the organism, simulation time was longer than 96 h in these scenarios.b Scenarios with p-values less than 0.05.

Table 6 .
Means of the criteria for model selection.

Table 6 .
Means of the criteria for model selection.

Table 6 .
Means of the criteria for model selection.

Table 7 .
Results of NLR and AUC estimation of representative drugs.

Table 7 .
Results of NLR and AUC estimation of representative drugs.

Table 8 .
Number of parameters required to characterize a plasma profile.